/*
** ShpTETR10.c  - this file contains methods for the 10-noded quadratic
**                tetrahedron solid finite element shape.
**
**    The adopted numbering scheeme is shown below:
**
**      The nodes are numbered according to the following rule:
**      1.  first, counterclockwise, the nodes at the face with 
**          natural coordinate L4 = 0,
**      2.  then, also counterclockwise, the nodes at the middle 
**          of the edges that go from the node with L4 = 1 to the 
**          the face with L4 = 0,
**      3.  and finally the node at the vertice with L4 = 1.
**
**                            + 10  
**                           /|\
**                          / | \
**                         /  |  \
**                        /   |   \
**                       /    |    \
**                     7/     |     \8
**                     /      |      \
**                    /       |9      \
**                   /        |        \
**                  /         |         \
**                 /          |          \
**               1+ . . . . . | .2. . . . +3
**                 \_         |         _/
**                   \_       |       _/
**                     \_     |     _/
**                      6\_   |   _/4
**                         \_ | _/
**                           \|/
**                            + 5                   
**
**
**           Internal numbering:
**
**
**                            9  (L4=1)
**                           /|\
**                          / | \
**                         /  |  \
**                        /   |   \
**                       /    |    \
**                      6     |     7
**                     /      |      \
**                    /       |       \
**                   /        8        \
**                  /         |         \
**                 /          |          \
**        (L1=1)  0 . . . . . | 1 . . . . 2  (L2=1)
**                 \_         |         _/
**                   \_       |       _/
**                     5_     |     _3
**                       \_   |   _/
**                         \_ | _/
**                           \|/
**                            4  (L3=1)
**
**
** --------------------------------------------------------------------------
**
** Entry Points:
**
** --------------------------------------------------------------------------
**
** void InitTETR10( void )
**
**      Initializes sub-class TETR10.
**
** --------------------------------------------------------------------------
**
** static Telement_3d *NewTETR10( void )
**
**      This method allocates memory for a given TETR10 element shape.
**
** --------------------------------------------------------------------------
**
** static void ReadReseqTETR10 ( int *reseq )
**
**      reseq    -  node numbering resequecing vector          (  in )
**
**      Reads node numbering resequecing vector TETR10 finite
**      element shape class.
**
** --------------------------------------------------------------------------
**
** static void ReseqTETR10 ( int *reseq )
**
**      reseq    -  node numbering resequecing vector          (  in )
**
**      Specifies node numbering resequecing for TETR10 finite
**      element shape class.
**
** --------------------------------------------------------------------------
**
** static void GetReseqTETR10 ( int *reseq )
**
**      reseq    -  node numbering resequecing vector          ( out )
**
**      Gets resequence vector for TETR10 finite element shape node
**      numbering.
**
** --------------------------------------------------------------------------
**
** static void ConnectTETR10 ( Telement_3d *shp,
**                             int    *connect )
**
**      shp      -  given element shape                        (  in )
**      connect  -  connectivity of the given element shape    ( out )
**
**      Returns connectivity of a given TETR10 finite element shape.
**
** --------------------------------------------------------------------------
**
** static void DerivMapTETR10rst ( Telement_3d    *shp,
**                                 NatCoord  *gauss_rst, 
**                                 DerivNat  deriv_rst[10] )
**
**      shp           -   given element shape  (not used)      (  in )
**      gauss_rst     -   natural coords.                      (  in )
**      deriv_rst     -   partial derivative wrt r,s,t         ( out )
**
**      This routine  returns the  "map"  functions  derivatives  with
**      respect to  the natural coordinates for a 10-noded tetrahedron
**      solid finite element shape. The "map" function derivatives are
**      evaluated at the given natural coordinates.
**      Numbering scheme is same as that given below for MapFuncTETR10.
**
** --------------------------------------------------------------------------
**
** static int  DimTETR10 ( Telement_3d *shp )
**
**      shp          -  given element shape  (not used)        (  in )
**
**      Returns dimension number of TETR10 shape, ie. 3.
**
**
** --------------------------------------------------------------------------
**
** static void GaussTypeTETR10 ( Telement_3d    *shp,
**                               GaussType *gauss_type )
**
**      shp        -  given element shape (not used)             (  in )
**      GaussType  -  returned Gauss integration type            ( out )
**
**      Returns the Gauss integration type which is consistent with 
**      this element shape, i.e., TETR_QUADRATURE.
**
** --------------------------------------------------------------------------
**
** static int  GetFirstEdgeTETR10 ( Telement_3d    *shp,
**                                  int       *num_verts,
**                                  int       connect[] )
**
**      shp        -  given element shape                        (  in )
**      num_verts  -  number edge verts                          ( out )
**      connect    -  edge connectivity                          ( out )
**
**      Gets no. of verts and connectivity of first edge of given shape.
**      The returned edge shape type is LINE3.
**
** --------------------------------------------------------------------------
**
** static int  GetNextEdgeTETR10 ( Telement_3d    *shp,
**                                 int       *num_verts,
**                                 int       connect[] )
**
**      shp        -  given element shape                        (  in )
**      num_verts  -  number edge verts                          ( out )
**      connect    -  edge connectivity                          ( out )
**
**      Gets no. of verts and connectivity of next edge of given shape.  
**      If there is no other edge returns a false (0) status.
**      The returned edge shape type is LINE3.
**
** --------------------------------------------------------------------------
**
** static int  GetFirstFaceTETR10 ( Telement_3d    *shp,
**                                  int       *num_verts,
**                                  int       connect[] )
**
**      shp        -  given element shape                        (  in )
**      num_verts  -  number face verts                          ( out )
**      connect    -  face connectivity                          ( out )
**
**      Gets no. of verts and connectivity of first face of given shape.
**      The returned face shape type is TRIA6.
**
** --------------------------------------------------------------------------
**
** static int  GetNextFaceTETR10 ( Telement_3d    *shp,
**                                 int       *num_verts,
**                                 int       connect[] )
**
**      shp        -  given element shape                        (  in )
**      num_verts  -  number face verts                          ( out )
**      connect    -  face connectivity                          ( out )
**
**      Gets no. of verts and connectivity of next face of given shape.  
**      If there is no other face returns a false (0) status.
**      The returned face shape type is TRIA6.
**
** --------------------------------------------------------------------------
**
** static int  GetShpEdgeTETR10  ( Telement_3d    *shp,
**                                 int       corners[],
**                                 ShapeType *type,
**                                 int       *num_nodes,
**                                 int       connect[] )
**
**      shp        -  given element shape                        (  in )
**      corners    -  corners of edge                            (  in )
**      type       -  returned edge shape type                   ( out )
**      num_nodes  -  number edge nodes                          ( out )
**      connect    -  edge connectivity                          ( out )
**
**      Gets shape type, no. of nodes, and connectivity of an edge of 
**      given shape for a given pair of corner nodes.
**      The returned edge shape type is LINE3.
**
** --------------------------------------------------------------------------
**
** static int  GetShpFaceTETR10  ( Telement_3d    *shp,
**                                 int       corners[],
**                                 ShapeType *type,
**                                 int       *num_nodes,
**                                 int       connect[] )
**
**      shp        -  given element shape                        (  in )
**      corners    -  corners of face                            (  in )
**      type       -  returned face shape type                   ( out )
**      num_nodes  -  number face nodes                          ( out )
**      connect    -  face connectivity                          ( out )
**
**      Gets shape type, no. of nodes, and connectivity of a face of 
**      given shape for a given triple of corner nodes.
**      The returned face shape type is TRIA6.
**
** --------------------------------------------------------------------------
**
** static void GetShpDomainTETR10  ( Telement_3d    *shp,
**                                   ShapeType *type,
**                                   int       *num_nodes,
**                                   int       connect[] )
**
**      shp        -  given element shape                        (  in )
**      type       -  returned domain shape type                 ( out )
**      num_nodes  -  number domain nodes                        ( out )
**      connect    -  domain connectivity                        ( out )
**
**      Gets shape type, no. of nodes, and connectivity of given shape.
**      The returned shape type is TETR10 and the number of domain nodes
**      is 10.
**
** --------------------------------------------------------------------------
**
** static void MakeTETR10  ( Telement_3d *shp,
**                           int    *connect )
**
**      shp      -  given element shape                        (  in )
**      connect  -  connectivity of the given element          (  in )
**
**      Makes a TETR10 element shape of given connectivity.
**
** --------------------------------------------------------------------------
**
** static void MakeDefaultTETR10  ( Telement_3d *shp,
**                                  int    *connect )
**
**      shp      -  given element shape                        (  in )
**      connect  -  connectivity of the given element          (  in )
**
**      Makes a TETR10 element shape of given connectivity.
**      The connectivity is assumed in the default class sequence.
**
** --------------------------------------------------------------------------
**
** static void MapFuncTETR10  ( Telement_3d    *shp,
**                              NatCoord  *nat_coord, 
**                              ShapeVal  *shape )
**
**      shp          -  given element shape (not used)         (  in )
**      nat_coord    -  natural coords.                        (  in )
**      shape        -  shape funcs. evaluated at given pt.    ( out )
**
**      This routine returns the "map" functions for a  10-noded solid
**      tetrahedron  element shape. The "map" functions are  evaluated
**      at the given natural coordinates. 
**
** --------------------------------------------------------------------------
**
** static void NodalNatCoordTETR10 ( Telement_3d    *shp,
**                                   NatCoord  *coord )
**
**      shp         -   given element shape                    (  in )
**      coords      -   natural coordinates of its nodes       ( out )
**
**      Returns shape natural nodal coordinates.
**
** --------------------------------------------------------------------------
**
** static int  NumNodesTETR10 ( Telement_3d *shp )
**
**      shp          -  given element shape  (not used)        (  in )
**
**      Returns number of nodes in a 10-noded tetrahedron shape, ie. 10.
**
** --------------------------------------------------------------------------
**
** static int  ReadTETR10 ( Telement_3d *shp, int *conn )
**
**      shp           -  Given element shape                   (  in )
**      conn          -  Connectivity of element               (  in )
**      
**      This method reads nodal connectivity  data of a given "TETR10" 
**      shape element.
**
** --------------------------------------------------------------------------
**
** version: 1-000
**
** history: Created:   30-Jul-96   Luiz Fernando Martha
**              Stolen from Femoop and modified.
**
**          Modified:  19-Oct-98   Joaquim Bento Cavalcante Neto
**              Modified some methods to agree with other shape modules.
**
** Modified:  Novem/2000 - Antonio Miranda (amiranda)
**      Changed the input of read method (FILE -> int *conn)
**
** --------------------------------------------------------------------------
*/

/* Global variables and symbols           */

#define _SHPTETR10_C

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#include "shp.h"


/*  Local variables and symbols  */

/*
** --------------------------------------------------------------------------
** Local data structures:
**
*/

typedef struct _datatetr10 {
   int   node[10];    /* element node indices */
} DataTETR10;

/* Resequence node numbering vector:
*/
static int  rs[10] = {  0,  1,  2,  3,  4,  5,  6,  7,  8,  9  };


/* Use quadrature 4 gauss point */
static int num_gauss_point = 4;
static NatCoord gauss_rst[4] = {
     { 0.5854101966249685, 0.1381966011250150, 0.1381966011250150 },
     { 0.1381966011250150, 0.5854101966249685, 0.1381966011250150 },
     { 0.1381966011250150, 0.1381966011250150, 0.5854101966249685 },
     { 0.1381966011250150, 0.1381966011250150, 0.1381966011250150 }
};
static GaussWgt gauss_wgt[4] = { 1.0/24.0, 1.0/24.0, 1.0/24.0, 1.0/24.0 };

static ShapeVal shape_rst[10][10];
static DerivNat deriv_rst[10][10];


/* Edge and Face lists:
 */
#define NEDGES  6
static int curr_edge = 0;;
static int el[NEDGES][3] =  { {   0,  1,  2  },
                  {   2,  3,  4  },
            {   4,  8,  9  },
            {   0,  5,  4  },
            {   0,  6,  9  },
            {   2,  7,  9  } };

#define NFACES  4
static int curr_face = 0;;
static int fl[NFACES][6] =  { {   0,  1,  2,  3,  4,  5  },
                              {   0,  5,  4,  8,  9,  6  },
            {   0,  6,  9,  7,  2,  1, },
            {   4,  3,  2,  7,  9,  8  } };

/* Local entry points: */

static Telement_3d *NewTETR10   ( void );
static void ReadReseqTETR10     ( int * );
static void ReseqTETR10         ( int * );
static void GetReseqTETR10      ( int * );
static void ConnectTETR10       ( Telement_3d *, int * );
static void DerivMapTETR10rst   ( Telement_3d *, NatCoord *, DerivNat[]  );
static int  DimTETR10           ( Telement_3d * );
static void MapFuncTETR10_i     ( Telement_3d *, int, ShapeVal **);

static int    NumGaussPoint           (void);
static void   DerivMapTETR10_i        ( Telement_3d *, int, DerivNat **);
static double GaussWgtTETR10_i        ( Telement_3d *shp, int i);
static void   GaussNatCoordTETR10_i   ( Telement_3d *shp, int i, NatCoord *);

static double NewSizeErrorTETR10      ( Telement_3d *, Tnode_data *, double);

static int  GetFirstEdgeTETR10  ( Telement_3d *, int *, int [] );
static int  GetNextEdgeTETR10   ( Telement_3d *, int *, int [] );
static int  GetFirstFaceTETR10  ( Telement_3d *, int *, int [] );
static int  GetNextFaceTETR10   ( Telement_3d *, int *, int [] );
static int  GetShpEdgeTETR10    ( Telement_3d *, int [], AuxShapeType *, int *, int [] );
static int  GetShpFaceTETR10    ( Telement_3d *, int [], AuxShapeType *, int *, int [] );
static void GetShpDomainTETR10  ( Telement_3d *, ShapeType *, int *, int[] );
static void MakeTETR10          ( Telement_3d *, int * );
static void MakeDefaultTETR10   ( Telement_3d *, int * );
static void MapFuncTETR10       ( Telement_3d *, NatCoord *, ShapeVal * );
static void NodalNatCoordTETR10 ( Telement_3d *, NatCoord * );
static int  NumNodesTETR10      ( Telement_3d * );
static int  ReadTETR10          ( Telement_3d *, int *conn );

/* Public Entry Points: */


/* ============================= InitTETR10 =============================== */

void InitTETR10( void )
{
   shape_class[NF_TETR10].New                   =    NewTETR10;
   shape_class[NF_TETR10].readresequence        =    ReadReseqTETR10;
   shape_class[NF_TETR10].resequence            =    ReseqTETR10;
   shape_class[NF_TETR10].getresequence         =    GetReseqTETR10;
   shape_class[NF_TETR10].connectivity          =    ConnectTETR10;
   shape_class[NF_TETR10].shapefunc             =    MapFuncTETR10;
   shape_class[NF_TETR10].shapefunc_i           =    MapFuncTETR10_i;

   shape_class[NF_TETR10].derivmaprst           =    DerivMapTETR10rst;
   shape_class[NF_TETR10].derivshprst           =    DerivMapTETR10rst;
   shape_class[NF_TETR10].derivmap_i            =    DerivMapTETR10_i;
   shape_class[NF_TETR10].derivshp_i            =    DerivMapTETR10_i;

   shape_class[NF_TETR10].derivxyz              =    shp_derivXYZ;
   shape_class[NF_TETR10].numgausspoint         =    NumGaussPoint;
   shape_class[NF_TETR10].bmatrix               =    shp_BMatrix;
   shape_class[NF_TETR10].gaussw_i              =    GaussWgtTETR10_i;
   shape_class[NF_TETR10].gaussnatcoord_i       =    GaussNatCoordTETR10_i;

   shape_class[NF_TETR10].dim                   =    DimTETR10;
   shape_class[NF_TETR10].free                  =    shp_freeshape;
   shape_class[NF_TETR10].gausstonode           =    shp_gausstonode;
   shape_class[NF_TETR10].jacobian              =    shp_jacobian;

   shape_class[NF_TETR10].getfirstedge          =    GetFirstEdgeTETR10;
   shape_class[NF_TETR10].getnextedge           =    GetNextEdgeTETR10;
   shape_class[NF_TETR10].getfirstface          =    GetFirstFaceTETR10;
   shape_class[NF_TETR10].getnextface           =    GetNextFaceTETR10;
   shape_class[NF_TETR10].getshpedge            =    GetShpEdgeTETR10;
   shape_class[NF_TETR10].getshpface            =    GetShpFaceTETR10;
   shape_class[NF_TETR10].getshpdomain          =    GetShpDomainTETR10;
   shape_class[NF_TETR10].getverts              =    ConnectTETR10;
   shape_class[NF_TETR10].isinterface           =    shp_isinterface;
   shape_class[NF_TETR10].localsys              =    shp_localsyssolid;
   shape_class[NF_TETR10].make                  =    MakeTETR10;
   shape_class[NF_TETR10].makedefault           =    MakeDefaultTETR10;
   shape_class[NF_TETR10].mapfunc               =    MapFuncTETR10;
   shape_class[NF_TETR10].nodalnatcoord         =    NodalNatCoordTETR10;
   shape_class[NF_TETR10].nodetovertex          =    shp_nodetovertex;
   shape_class[NF_TETR10].nummapnodes           =    NumNodesTETR10;
   shape_class[NF_TETR10].numshpnodes           =    NumNodesTETR10;
   shape_class[NF_TETR10].numverts              =    NumNodesTETR10;
   shape_class[NF_TETR10].print                 =    shp_printshape;
   shape_class[NF_TETR10].read                  =    ReadTETR10;

   shape_class[NF_TETR10].sizeerror             =    NewSizeErrorTETR10;

   int i;
   for (i = 0; i < num_gauss_point; i++)
   {
     MapFuncTETR10 (NULL, &gauss_rst[i], shape_rst[i]);
     DerivMapTETR10rst ( NULL, &gauss_rst[i], deriv_rst[i]);
   }   
}

/* ============================== NewTETR10 ================================= */

static Telement_3d *NewTETR10 ( void )
{
   Telement_3d       *shp;

   shp = (Telement_3d *)calloc( 1, sizeof(Telement_3d) );

   shp->type = NF_TETR10;

   return( shp );
}


/* ============================ ReadReseqTETR10 ============================= */

static void ReadReseqTETR10 ( int *reseq )
{
   int i;

   for (i=0;i<10;i++)
    rs[i] = reseq[i] - 1;
}


/* ============================== ReseqTETR10 =============================== */

static void ReseqTETR10 ( int *reseq )
{
   int i;

   for( i=0; i<10; i++ )
    rs[i] = reseq[i];
}


/* ============================= GetReseqTETR10 ============================= */

static void GetReseqTETR10 ( int *reseq )
{
   int i;

   for( i=0; i<10; i++ )
     reseq[i] = rs[i];
}

/* =============================== ConnectTETR10 ============================ */

static void ConnectTETR10 ( Telement_3d *shp,
                            int    *connect )
{
   connect[0] = shp->inc[0];
   connect[1] = shp->inc[1];
   connect[2] = shp->inc[2];
   connect[3] = shp->inc[3];
   connect[4] = shp->inc[4];
   connect[5] = shp->inc[5];
   connect[6] = shp->inc[6];
   connect[7] = shp->inc[7];
   connect[8] = shp->inc[8];
   connect[9] = shp->inc[9];
}


/* ========================== DerivMapTETR10rst ============================= */

static void DerivMapTETR10rst ( Telement_3d   *shp,
                                NatCoord *gauss_rst,
                                DerivNat deriv_rst[10] )
{
   double  L1, L2, L3, L4;

   ignore (shp);

   L1 = gauss_rst->r;
   L2 = gauss_rst->s;
   L3 = gauss_rst->t;
   L4 = 1.0 - L1 - L2 - L3;


   deriv_rst[0].r =  4.0 * L1 - 1.0;
   deriv_rst[1].r =  4.0 * L2;
   deriv_rst[2].r =  0.0;
   deriv_rst[3].r =  0.0;
   deriv_rst[4].r =  0.0;
   deriv_rst[5].r =  4.0 * L3;
   deriv_rst[6].r =  4.0 * (L4-L1);
   deriv_rst[7].r = -4.0 * L2;
   deriv_rst[8].r = -4.0 * L3;
   deriv_rst[9].r =  1.0 - 4.0 * L4;
   

   deriv_rst[0].s =  0.0;
   deriv_rst[1].s =  4.0 * L1;
   deriv_rst[2].s =  4.0 * L2 - 1.0;
   deriv_rst[3].s =  4.0 * L3;
   deriv_rst[4].s =  0.0;
   deriv_rst[5].s =  0.0;
   deriv_rst[6].s = -4.0 * L1;
   deriv_rst[7].s =  4.0 * (L4-L2);
   deriv_rst[8].s = -4.0 * L3;
   deriv_rst[9].s =  1.0 - 4.0 * L4;
   

   deriv_rst[0].t =  0.0;
   deriv_rst[1].t =  0.0;
   deriv_rst[2].t =  0.0;
   deriv_rst[3].t =  4.0 * L2;
   deriv_rst[4].t =  4.0 * L3 - 1.0;
   deriv_rst[5].t =  4.0 * L1;
   deriv_rst[6].t = -4.0 * L1;
   deriv_rst[7].t = -4.0 * L2;
   deriv_rst[8].t =  4.0 * (L4-L3);
   deriv_rst[9].t =  1.0 - 4.0 * L4;
}

/* ============================== NumGaussPoint =========================== */
static int NumGaussPoint (void)
{
  return (num_gauss_point);
}

/* ============================== DerivMapTETRA10_i =========================== */
static void DerivMapTETR10_i  ( Telement_3d *shp, int i, DerivNat **ptr)
{
  if (i < 0 || i > 9)
  {
    ptr = NULL;
    return;
  }
  *ptr = deriv_rst[i];
}


/* ============================== GaussWgtTETRA10_i =========================== */

static double GaussWgtTETR10_i  ( Telement_3d *shp, int i)
{
  if (i < 0 || i > 4)
    return 0.0;

  return (gauss_wgt[i]);
}


/* ============================== GaussNatCoordTETR10_i =========================== */

static void GaussNatCoordTETR10_i ( Telement_3d *shp, int i, NatCoord *natCoord)
{
  if (i < 0 || i > 4)
    return;
  natCoord->r = gauss_rst[i].r;
  natCoord->s = gauss_rst[i].s;
  natCoord->t = gauss_rst[i].t;
}


/* =============================== DimTETR10 ================================ */

static int DimTETR10( Telement_3d *shp )
{
   ignore (shp);

   return ( 3 );
}


/* ============================= GetFirstEdgeTETR10 ========================= */

static int  GetFirstEdgeTETR10  ( Telement_3d    *shp,
           int       *num_verts,
                 int       connect[] )
{
   *num_verts = 3;
   curr_edge = 0;

   connect[0] = shp->inc[el[curr_edge][0]];
   connect[1] = shp->inc[el[curr_edge][1]];
   connect[2] = shp->inc[el[curr_edge][2]];
   return( 1 );
}


/* ============================= GetNextEdgeTETR10 ========================== */

static int  GetNextEdgeTETR10 ( Telement_3d    *shp,
              int       *num_verts,
              int       connect[] )
{
   *num_verts = 3;
   curr_edge++;

   if( curr_edge >= NEDGES )
    return( 0 );

   connect[0] = shp->inc[el[curr_edge][0]];
   connect[1] = shp->inc[el[curr_edge][1]];
   connect[2] = shp->inc[el[curr_edge][2]];
   return( 1 );
}


/* ============================= GetFirstFaceTETR10 ======================== */

static int  GetFirstFaceTETR10 ( Telement_3d    *shp,
         int       *num_verts,
         int       connect[] )
{
   *num_verts = 6;
   curr_face = 0;

   connect[0] = shp->inc[fl[curr_face][0]];
   connect[1] = shp->inc[fl[curr_face][1]];
   connect[2] = shp->inc[fl[curr_face][2]];
   connect[3] = shp->inc[fl[curr_face][3]];
   connect[4] = shp->inc[fl[curr_face][4]];
   connect[5] = shp->inc[fl[curr_face][5]];
   return( 1 );
}


/* ============================= GetNextFaceTETR10 =========================== */

static int  GetNextFaceTETR10  ( Telement_3d    *shp,
        int       *num_verts,
        int       connect[] )
{
   *num_verts = 6;
   curr_face++;

   if( curr_face >= NFACES )
    return( 0 );

   connect[0] = shp->inc[fl[curr_face][0]];
   connect[1] = shp->inc[fl[curr_face][1]];
   connect[2] = shp->inc[fl[curr_face][2]];
   connect[3] = shp->inc[fl[curr_face][3]];
   connect[4] = shp->inc[fl[curr_face][4]];
   connect[5] = shp->inc[fl[curr_face][5]];
   return( 1 );
}
/* =============================== GetShpEdgeTETR10 ========================= */

static int  GetShpEdgeTETR10 ( Telement_3d    *shp,
                               int       corners[],
                               AuxShapeType *type,
                               int       *num_nodes,
                               int       connect[] )
{
   int          status = 0;

   *type = LINE3;
   *num_nodes = 3;

   if( corners[0] == shp->inc[0] )
   {
    connect[0] = shp->inc[0];
    if( corners[1] == shp->inc[9] )
    {
     connect[1] = shp->inc[6];
     connect[2] = shp->inc[9];
     status = 1;
    }
    else if( corners[1] == shp->inc[2] )
    {
     connect[1] = shp->inc[1];
     connect[2] = shp->inc[2];
     status = 1;
    }
    else if( corners[1] == shp->inc[4] )
    {
     connect[1] = shp->inc[5];
     connect[2] = shp->inc[4];
     status = 1;
    }
   }
   else if( corners[0] == shp->inc[2] )
   {
    connect[0] = shp->inc[2];
    if( corners[1] == shp->inc[0] )
    {
     connect[1] = shp->inc[1];
     connect[2] = shp->inc[0];
     status = 1;
    }
    else if( corners[1] == shp->inc[4] )
    {
     connect[1] = shp->inc[3];
     connect[2] = shp->inc[4];
     status = 1;
    }
    else if( corners[1] == shp->inc[9] )
    {
     connect[1] = shp->inc[7];
     connect[2] = shp->inc[9];
     status = 1;
    }
   }
   else if( corners[0] == shp->inc[4] )
   {
    connect[0] = shp->inc[4];
    if( corners[1] == shp->inc[2] )
    {
     connect[1] = shp->inc[3];
     connect[2] = shp->inc[2];
     status = 1;
    }
    else if( corners[1] == shp->inc[0] )
    {
     connect[1] = shp->inc[5];
     connect[2] = shp->inc[0];
     status = 1;
    }
    else if( corners[1] == shp->inc[9] )
    {
     connect[1] = shp->inc[8];
     connect[2] = shp->inc[9];
     status = 1;
    }
   }
   else if( corners[0] == shp->inc[9] )
   {
    connect[0] = shp->inc[9];
    if( corners[1] == shp->inc[4] )
    {
     connect[1] = shp->inc[8];
     connect[2] = shp->inc[4];
     status = 1;
    }
    else if( corners[1] == shp->inc[0] )
    {
     connect[1] = shp->inc[6];
     connect[2] = shp->inc[0];
     status = 1;
    }
    else if( corners[1] == shp->inc[2] )
    {
     connect[1] = shp->inc[7];
     connect[2] = shp->inc[2];
     status = 1;
    }
   }
   
   return( status );
}


/* =============================== GetShpFaceTETR10 ========================= */

static int  GetShpFaceTETR10 ( Telement_3d    *shp,
                               int       corners[],
                               AuxShapeType *type,
                               int       *num_nodes,
                               int       connect[] )
{
   int          status = 0;

   *type = TRIA6;
   *num_nodes = 6;

   if( corners[0] == shp->inc[0] )
   {
    connect[0] = shp->inc[0];
    if( (corners[1] == shp->inc[2]) && (corners[2] == shp->inc[4]) )
    {
     connect[1] = shp->inc[1];
     connect[2] = shp->inc[2];
     connect[3] = shp->inc[3];
     connect[4] = shp->inc[4];
     connect[5] = shp->inc[5];
     status = 1;
    }
    else if( (corners[1] == shp->inc[4]) && (corners[2] == shp->inc[9]) )
    {
     connect[1] = shp->inc[5];
     connect[2] = shp->inc[4];
     connect[3] = shp->inc[8];
     connect[4] = shp->inc[9];
     connect[5] = shp->inc[6];
     status = 1;
    }
    else if( (corners[1] == shp->inc[9]) && (corners[2] == shp->inc[2]) )
    {
     connect[1] = shp->inc[6];
     connect[2] = shp->inc[9];
     connect[3] = shp->inc[7];
     connect[4] = shp->inc[2];
     connect[5] = shp->inc[1];
     status = 1;
    }
   }
   else if( corners[0] == shp->inc[2] )
   {
    connect[0] = shp->inc[2];
    if( (corners[1] == shp->inc[9]) && (corners[2] == shp->inc[4]) )
    {
     connect[1] = shp->inc[7];
     connect[2] = shp->inc[9];
     connect[3] = shp->inc[8];
     connect[4] = shp->inc[4];
     connect[5] = shp->inc[3];
     status = 1;
    }
    else if( (corners[1] == shp->inc[0]) && (corners[2] == shp->inc[9]) )
    {
     connect[1] = shp->inc[1];
     connect[2] = shp->inc[0];
     connect[3] = shp->inc[6];
     connect[4] = shp->inc[9];
     connect[5] = shp->inc[7];
     status = 1;
    }
    else if( (corners[1] == shp->inc[4]) && (corners[2] == shp->inc[0]) )
    {
     connect[1] = shp->inc[3];
     connect[2] = shp->inc[4];
     connect[3] = shp->inc[5];
     connect[4] = shp->inc[0];
     connect[5] = shp->inc[1];
     status = 1;
    }
   }
   else if( corners[0] == shp->inc[4] )
   {
    connect[0] = shp->inc[4];
    if( (corners[1] == shp->inc[0]) && (corners[2] == shp->inc[2]) )
    {
     connect[1] = shp->inc[5];
     connect[2] = shp->inc[0];
     connect[3] = shp->inc[1];
     connect[4] = shp->inc[2];
     connect[5] = shp->inc[3];
     status = 1;
    }
    else if( (corners[1] == shp->inc[2]) && (corners[2] == shp->inc[9]) )
    {
     connect[1] = shp->inc[3];
     connect[2] = shp->inc[2];
     connect[3] = shp->inc[7];
     connect[4] = shp->inc[9];
     connect[5] = shp->inc[8];
     status = 1;
    }
    else if( (corners[1] == shp->inc[9]) && (corners[2] == shp->inc[0]) )
    {
     connect[1] = shp->inc[8];
     connect[2] = shp->inc[9];
     connect[3] = shp->inc[6];
     connect[4] = shp->inc[0];
     connect[5] = shp->inc[5];
     status = 1;
    }
   }
   else if( corners[0] == shp->inc[9] )
   {
    connect[0] = shp->inc[9];
    if( (corners[1] == shp->inc[0]) && (corners[2] == shp->inc[4]) )
    {
     connect[1] = shp->inc[6];
     connect[2] = shp->inc[0];
     connect[3] = shp->inc[5];
     connect[4] = shp->inc[4];
     connect[5] = shp->inc[8];
     status = 1;
    }
    else if( (corners[1] == shp->inc[4]) && (corners[2] == shp->inc[2]) )
    {
     connect[1] = shp->inc[8];
     connect[2] = shp->inc[4];
     connect[3] = shp->inc[3];
     connect[4] = shp->inc[2];
     connect[5] = shp->inc[7];
     status = 1;
    }
    else if( (corners[1] == shp->inc[2]) && (corners[2] == shp->inc[0]) )
    {
     connect[1] = shp->inc[7];
     connect[2] = shp->inc[2];
     connect[3] = shp->inc[1];
     connect[4] = shp->inc[0];
     connect[5] = shp->inc[6];
     status = 1;
    }
   }

   return( status );
}


/* ============================== GetShpDomainTETR10 ======================== */

static void GetShpDomainTETR10 ( Telement_3d    *shp,
                                 ShapeType *type,
                                 int       *num_nodes,
                                 int       connect[] )
{
   *type = NF_TETR10;
   *num_nodes = 10;

   ConnectTETR10( shp, connect );
}


/* ================================ MakeTETR10 ============================== */

static void MakeTETR10 ( Telement_3d *shp,
                         int    *connect )
{
   shp->inc[rs[0]] = connect[0];
   shp->inc[rs[1]] = connect[1];
   shp->inc[rs[2]] = connect[2];
   shp->inc[rs[3]] = connect[3];
   shp->inc[rs[4]] = connect[4];
   shp->inc[rs[5]] = connect[5];
   shp->inc[rs[6]] = connect[6];
   shp->inc[rs[7]] = connect[7];
   shp->inc[rs[8]] = connect[8];
   shp->inc[rs[9]] = connect[9];
}

/* ================================ MakeDefaultTETR10 ======================= */

static void MakeDefaultTETR10 ( Telement_3d *shp,
                                int    *connect )
{
   shp->inc[0] = connect[0];
   shp->inc[1] = connect[1];
   shp->inc[2] = connect[2];
   shp->inc[3] = connect[3];
   shp->inc[4] = connect[4];
   shp->inc[5] = connect[5];
   shp->inc[6] = connect[6];
   shp->inc[7] = connect[7];
   shp->inc[8] = connect[8];
   shp->inc[9] = connect[9];
}


/* ============================= MapFuncTETR10 ============================== */

static void MapFuncTETR10 ( Telement_3d   *shp,
                            NatCoord *nat_coord, 
                            ShapeVal *shape )
{
   double  L1, L2, L3, L4;

   ignore (shp);

   L1 = nat_coord->r;
   L2 = nat_coord->s;
   L3 = nat_coord->t;
   L4 = 1.0 - L1 - L2 - L3;

 
   shape[0] = (2.0 * L1 - 1.0) * L1;
   shape[1] =  4.0 * L1 * L2;
   shape[2] = (2.0 * L2 - 1.0) * L2;
   shape[3] =  4.0 * L2 * L3;
   shape[4] = (2.0 * L3 - 1.0) * L3;
   shape[5] =  4.0 * L1 * L3;
   shape[6] =  4.0 * L1 * L4;
   shape[7] =  4.0 * L2 * L4;
   shape[8] =  4.0 * L3 * L4;
   shape[9] = (2.0 * L4 - 1.0) * L4;
}

/* ============================= MapFuncTETR10_i ============================== */

static void MapFuncTETR10_i (Telement_3d *shp, int i, ShapeVal **shape )
{
  if (i < 0 || i > num_gauss_point)
  {
    *shape = NULL;
    return;
  }
  *shape = shape_rst[i];
}

/* =========================== NodalNatCoordTETR10 ========================== */

static void NodalNatCoordTETR10 ( Telement_3d    *shp,
                                  NatCoord  *coord )
{
   ignore (shp);

   coord[0].r =  1.0; coord[0].s =  0.0; coord[0].t =  0.0;
   coord[1].r =  0.5; coord[1].s =  0.5; coord[1].t =  0.0;
   coord[2].r =  0.0; coord[2].s =  1.0; coord[2].t =  0.0;
   coord[3].r =  0.0; coord[3].s =  0.5; coord[3].t =  0.5;
   coord[4].r =  0.0; coord[4].s =  0.0; coord[4].t =  1.0;
   coord[5].r =  0.5; coord[5].s =  0.0; coord[5].t =  0.5;
   coord[6].r =  0.5; coord[6].s =  0.0; coord[6].t =  0.0;
   coord[7].r =  0.0; coord[7].s =  0.5; coord[7].t =  0.0;
   coord[8].r =  0.0; coord[8].s =  0.0; coord[8].t =  0.5;
   coord[9].r =  0.0; coord[9].s =  0.0; coord[9].t =  0.0;

}


/* ============================ NumNodesTETR10 ============================== */

static int NumNodesTETR10( Telement_3d *shp )
{
   ignore (shp);

   return ( 10 );
}


/* ============================= ReadTETR10 ================================= */

static int ReadTETR10( Telement_3d *shp, int *conn )
{
   int i;

   for (i=0;i<10;i++)
     shp->inc[rs[i]] = conn[i];

   return( TRUE );
}


/* =========================== NewSizeErrorTETR10 ============================= */

static double NewSizeErrorTETR10 (Telement_3d *shp, Tnode_data *node_coords, 
                                  double ratio_error)
{
  double  u[3], v[3], w[3];
  double  vol_real;
  double  hold;
  double  hnew;
  double  fact;
  int      i;
  Tnode_data coord[10];

  // obtain node coords
  for (i = 0; i < shp->n; ++i)
  {
    coord[i].x = node_coords[shp->inc[i]].x;
    coord[i].y = node_coords[shp->inc[i]].y;
    coord[i].z = node_coords[shp->inc[i]].z;
  }

  // fill the three component vectors that define  the tetraedron
  u[0] = coord[2].x - coord[0].x;
  u[1] = coord[2].y - coord[0].y;
  u[2] = coord[2].z - coord[0].z;
  v[0] = coord[4].x - coord[0].x;
  v[1] = coord[4].y - coord[0].y;
  v[2] = coord[4].z - coord[0].z;
  w[0] = coord[9].x - coord[0].x;
  w[1] = coord[9].y - coord[0].y;
  w[2] = coord[9].z - coord[0].z;

  // get the actual element characteristic size
  vol_real = (1.0/6.0)*fabs((u[0]*v[1]*w[2] + u[2]*v[0]*w[1] + u[1]*w[0]*v[2])-
                            (w[0]*v[1]*u[2] + v[0]*u[1]*w[2] + u[0]*v[2]*w[1]));

  hold = pow((12.0*vol_real)/(sqrt(2.0)),(1.0/3.0));

  // calculate the new element size
  fact = pow (ratio_error, 1.0/2.0 );

  if( fact <= 0.0 ) 
    hnew = hold;
  else
    hnew = hold / fact;

  // return the new element size based in a error estimative
  return (hnew);

}

/* =========================================================== End of File */
